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ABSTRACT 

Here we present a new approach for constraining luminous blazars, incorporating fully time- 
dependent and self-consistent modeling of bright 7 -ray flares of PKS 1510-089 resolved with Fermi- 
LAT, in the framework of the internal shock scenario. The results of our modeling imply the location 
of the 7 -ray flaring zone outside of the broad-line region, namely around ~ 0.3 pc from the core for 
a free-expanding jet with the opening angle T 0 jet — 1 (where P is the jet bulk Lorentz factor), up 
to ~ 3 pc for a collimated outflow with r djet — 0 . 1 . Moreover, under the r 9jet — 1 condition, our 
modeling indicates the maximum efficiency of the jet production during the flares, with the total 
jet energy flux strongly dominated by protons and exceeding the available accretion power in the 
source. This is in contrast to the quiescence states of the blazar, characterized by lower jet kinetic 
power and an approximate energy equipartition between different plasma constituents. We demostrate 
how strictly simultaneous observations of flaring PKS 1510-089 at optical, X-ray, and GeV photon 
energies on hourly timescales, augmented by extensive simulations as presented in this paper, may 
help to impose further precise constraints on the magnetization and opening angle of the emitting 
region. Our detailed modeling implies in addition that a non-uniformity of the Doppler factor across 
the jet, caused by the radial expansion of the outflow, may lead to a pronounced time distortion in 
the observed 7 -ray light curves, resulting in particular in asymmetric flux profiles with substantially 
extended decay phases. 

Subject headings: acceleration of particles — radiation mechanisms: non-thermal — galaxies: active 
— galaxies: jets — quasars: individual (PKS 1510—089) — gamma rays: galaxies 


1. INTRODUCTION 

Blazars constitute a population of radio-loud Active 
Galactic Nuclei (AGN) with relativistic jets directed to¬ 
ward the Earth. The observed spectrum of a blazar is 
dominated by the jet component as a result of a strong 
relativistic beaming involved, and consists of the two 
broad humps in the v — vFi, representation (hereafter 
the ‘spectral energy distribution’, SED). The low energy 
hump, peaked in the infrared-to-X-ray frequency range, 
is due to the synchrotron emission of ultra-relativistic 
jet electrons, while the high energy component, peaking 
in the 7 -ray band, originates from the inverse Gompton 
radiation of relativistic electrons interacting with ambi¬ 
ent seed photons produced either externally to the jet 
(predominantly via a re-procession of the accretion disk 
emission by the circumnuclear gas and dust), or within 
the outflow by the synchrotron process. 

Blazars can be classified as either BL Lacertae objects 
(BL Lacs) or flat-spectrum radio quasars (FSRQs), based 
on the properties of the optical line emission. In partic¬ 
ular, the latter class of sources is characterized by the 
presence of prominent broad and narrow emission lines 
in their spectra. FSRQs constitute the most luminous 
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population of blazars, with super massive black holes 
(SMBHs) accreting at high rates (over 1% of the Ed¬ 
dington limit; e.g., Sbarrato et al. 2014). The available 
huge accretion mass flux is converted very efficiently to 
the jet energy flux within the ergospheres of rapidly spin¬ 
ning SMBHs at the expense of the black hole rotational 
energy (Blandford & Znajek 1977; Sikora et al. 2007; 
Tchekhovskoy et al. 2011). In this context, high-energy 
observations are indispensable for imposing meaningful 
constraints on the jet energetics, since the bulk of the 
radiative energy of luminous blazars is released in the 
MeV/GeV range (see Ghisellini et al. 2014). And in¬ 
deed, it was repeatedly shown that the amount of the jet 
power dissipated radiatively during the flares of FSRQs 
in 7 -rays is often comparable to the corresponding disk 
luminosities (Tanaka et al. 2011; Saito et al. 2013). 

Blazars display complex variability patterns on var¬ 
ious timescales across the entire electromagnetic spec¬ 
trum, with high-amplitude flux changes pronounced 
most clearly at high frequencies. The all-sky 7 -ray survey 
with the Large Area Telescope onboard the Fermi satel¬ 
lite (hereafter Fermi-LAT) unveiled in particular very 
rapid flares of FSRQs in the GeV range, with the ob¬ 
served flux doubling timescales of the order of hours, 
down to even sub-hour domain in a few cases (see the 
analysis of the LAT data for PKS 1510-089, 3C273, and 
4C 21.35 by Foschini et al. 2011, 2013; Saito et al. 2013; 
Brown 2013; Rani et al. 2013). Such a dramatic 7 -ray 
variability indicates a very efficient energization of the jet 
electrons taking place in compact regions of the enhanced 
energy dissipation, followed by the rapid radiative cool¬ 
ing of the accelerated particles. 

The exact site and structure of those energy dissipa- 
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tion regions is, however, still controversial, in spite of 
the extensive broad-band observations carried out over 
the last two decades. This reflects in a wide range of 
the estimated position of the blazar emission zone along 
the jet, even for the same object, based on different sets 
of observables interpreted in the frameworks of differ¬ 
ent jet emission models. For example, as for PKS 1510-- 
089, the combined radio, optical polarization, and 7 -ray 
monitoring suggested the dominant emitting region to 
be located as far as ~ 10 pc from the central SMBH (the 
“far-dissipation zone” scenario; e.g., Marscher et al. 2010; 
Orient! et al. 2013). Meanwhile, another studies based 
on the SED modeling including the UV, X-ray, and 7 - 
ray data, placed the blazar zone in the source at sub-pc 
distances from the central engine (the “near-dissipation 
zone” scenario; e.g., Kataoka et al. 2008; Ghisellini et 
al. 2010). Complex multi-zone/multi-component models 
were also proposed in order to explain broad-band vari¬ 
ability and spectral properties of PKS 1510-089 (Nale- 
wajko et al. 2012a; Aleksic et al. 2014). 

In principle, timing characteristics of high-amplitude 
blazar flares can offer some clues on the location of the 
active emission zone. First, due to the causality re¬ 
quirement, the linear size of the emission region may 
be constrained as r < c6Td/{l + z), for the observed 
flux doubling timescale Td and the Doppler factor of the 
emitting plasma S. Assuming further a conical, uniform, 
and free-expanding jet, whose (small) half-opening angle 
0 jet — r/R is approximately the inverse of the jet bulk 
Lorentz factor F, the location of the flaring zone from 
the central SMBH reads as i? < cF STd/{l + z). The flux 
doubling timescale of about a few hours, for example, of¬ 
ten inferred for FSRQs from the LAT data as mentioned 
above, gives therefore R < 10^® cm for the typically ex¬ 
pected in blazar sources F ~ ~ 10, or i? < 10^ Rs in 

the units of the Schwarzschild radius Rs for the black 
hole mass of the order of 10® M 0 . This would be then 
consistent with the “near-dissipation zone” scenario. 

On the other hand, the most recent detections of very 
high energy 7 -rays (photon energies > 100 GeV) from 
some FSRQs, including PKS 1510-089 (Abramowski et 
al. 2013; Aleksic et al. 2014; Barnacka et al. 2014), have 
challenged the above conclusion. That is because such 
7 -rays are expected to undergo an efficient annihilation 
in the interactions with lower-energy circumnuclear pho¬ 
ton fields, and in particular with the UV disk emission 
re-processed within the “broad line region” (BLR), if pro¬ 
duced at the scales below the characteristic radius of the 
BLR clouds (~ 0.1 pc for the typical disk luminosity of 
the order of lO'^^ergs”^; see, e.g., Moderski et al. 2003; 
Liu & Bai 2006; Tavecchio & Mazin 2009; Barnacka et 
al. 2014). The detection of very high energy 7 -rays from 
FSRQs seems therefore to require the dominant emission 
region to be located at larger distances from the core, at 
least in the framework of the one-zone emission models. 

An important issue, which is currently overlooked in 
the blazar modeling, but which, in fact, is crucial for 
constraining luminous blazars, is the exact 7 -ray spec¬ 
tral evolution during the rapid flux changes, along with 
the exact time profiles of the flares. Fermi-LAT observa¬ 
tions of the brightest FSRQs have revealed that spectral 
breaks appear around a few/several GeV during the en¬ 
hanced activity states, and that the low- and high-energy 


TABLE 1 

7-RAY FLARING EVENTS OF PKS 1510-089 MODELED IN THIS WORK 


Name 

MJD 

^>100 MeV 



( 1 ) 

( 2 ) 

(3) 

Flare 1 

55853.5-55854.5 

14.86 ± 0.89 

1.97 ± 0.04 

Flare 2 

55872-55874 

8.39 ± 0.44 

2.19 ± 0.04 


(1) Dates of the 7 -ray flux maxima in the daily-binned light curve; 

( 2 ) photon fluxes measured at the flux maxima in the units of 

[10“® phcm“^ averaged over the specified time intervals; (3) 

the corresponding photon indices. All the values are taken from 
Saito et al. (2013). 

spectral slopes below and above the break, respectively, 
often change significantly on the timescales of less than a 
day (Abdo et al. 2009, 2011; Tanaka et al. 2011; Stern & 
Poutanen 2011; Rani et al. 2013). This indicates that the 
7 -ray data integrated over longer periods of time, which 
are typically utilized in the blazar SED modeling, should 
be considered as a superposition of different spectra pro¬ 
duced at different stages of the jet evolution, and as such 
may hardly be used for any exact spectral diagnosis. Of 
course, in the overwhelming majority of cases long in¬ 
tegration LAT exposures times are unavoidable, due to 
a limited photon statistics. Still, the caution should be 
kept in mind. 

This situation has motivated us to attempt a detailed, 
time-dependent, and fully self-consistent modeling of the 
spectral evolution of bright 7 -ray flares in ESRQs, for 
which the time profiles can be resolved and the spectra 
can be constrained on the hourly timescales. In partic¬ 
ular, in this paper we present the extensive modeling of 
the extremely bright 7 -ray flares in PKS 1510-089, which 
were analyzed before by Saito et al. (2013, see also Brown 
2013 and Foschini et al. 2013). This allows us to put 
uniquely robust (though still mo del-dependent, to some 
extent) constraints on the structure and the location of 
the flaring emission zone in the source, as discussed fur¬ 
ther below. 

Throughout the paper we assume the AGDM cos¬ 
mology with Da = 0.73, Dm = 0.27, and Hq = 
71 kms“^ Mpc“^, so that the redshift of PKS 1510-089, 
z = 0.361, corresponds to the luminosity distance of 
c?L — 1.91 Gpc. 

2. SELECTION OF THE FERMI-LAT DATA 

The Fermi-LAT provides the most complete and sen¬ 
sitive coverage of the 7 -ray sky up to date. Since 2008 
it surveys the entire sky every three hours within the 
photon energy range 20 MeV — 300 GeV, so that each as- 
trophysical 7 -ray source — majority of which are blazars 
(Nolan et al. 2012) — is exposed for about 30 min during 
one scan. Due to the limited photon statistics, the selec¬ 
tion of the Fermi-LAT data for this work has to focus in¬ 
evitably on the flaring states of the brightest FSRQs, for 
which meaningful spectral analysis could be performed 
with the minimum LAT time resolution. We further re¬ 
strict our modeling to the well-resolved 7 -ray flares of 
PKS 1510-089 which occurred in 2011, as summarized 
in Table 1 following the analysis of Saito et al. (2013). 
These two flares constitute the best known examples of 
the prominent, isolated, and coherent events, unlike the 
majority of the observed blazar flux enhancements which 
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TABLE 2 


Model parameters oe the “eree-expanding jet” fit to the 7-ray flares of PKS 1510-089. 

Model parameter 

Flare 1 Flare 2 

reference 

Minimum electron Lorentz factor, 7min 


1 

Barnacka et al. (2014) 

Break electron Lorentz factor, Tbr 


900 


Maximum electron Lorentz factor, 7max 


1 X 10® 


Low-energy electron injection index, p 


1.2 


High-energy electron injection index, q 


3.4 

” 

Bulk Lorentz factor of the emitting shell, T 


22 

55 

Jet opening angle, 0jet 


2.6 deg 

55 

Jet viewing angle, ^obs 


2.6 deg 

55 

Jet magnetic field intensity at 10^® cm, Bq 


0.75 G 

55 

Characteristic scale of the BLR, -Rblr 


0.12 X 10^® cm 

55 

Central energy density of the BLR photon field, wblr 


0.06 ergcm”^ 

55 

Characteristic energy of the BLR photons, Ai^blr 


10 eV 

55 

Characteristic scale of the HDT, Rhdt 


1.94 X 10l®cm 

55 

Central energy density of the HDT photon field, whdt 


5 X 10“"* ergcm“® 

55 

Characteristic energy of the HDT photons, Az/hdt 


0.15eV 

55 

Normalization of the electron injection function, Ke 

1.6 X 

lO'i'^s-i 0.6 X lO'^’^s-i 

this work, § 4.1 

Distance where the injection starts, Rstart 

0.7 X 

10 l®cm 2.3x101® cm 

55 

Distance where the injection terminates, -Rstop 

0.9 X 

lOl® cm 3.4 X lOl® cm 

55 

Distance where the simulation stops, Rend 

2.3 X 

lOl® cm 6.9xl0l®cm 

55 


seem rather like a superposition of distinct (though pos¬ 
sibly related) but just unresolved sub-flaring/flickering. 

The question remains if the minimum (orbital) 3-hour 
binning of the LAT data is sufficiently short for resolving 
the PKS 1510-089 flares properly. We note in this con¬ 
text that the stochastic modeling of the LAT light curves 
for the brightest blazars by Sobolewska et al. (2014) re¬ 
vealed some hints for sub-hour variability timescales in 
only four sources, not including PKS 1510-089. These 
characteristic variability timescales regarded the features 
of the 7 -ray power spectral density functions, and not the 
flux doubling timescales studied by Saito et al. (2013) or 
Brown (2013). Taking therefore into account that the 
rising and decaying phases of the two flares selected for 
our analysis are very well defined, with no sub-structure 
appearing even when studied with the shortest binning 
(Saito et al. 2013), we conclude that one can indeed con¬ 
sider a single emitting component to be responsible for 
the observed flux evolution during the targeted activity 
epochs of the source, and hence that a meaningful time- 
dependent modeling can be performed in the framework 
of a given jet model. The issue of sub-orbital blazar vari¬ 
ability in the LAT data will be analyzed in detail in the 
forthcoming paper Saito et al. (2015, in prep.). 

3. THE MODEL 

We simulate the time evolution of the observed source 
spectrum during the selected flaring events utilizing the 
BLAZAR code, which was developed by Moderski et al. 
(2003, 2005) based on the widely anticipated scenario 
of internal shocks formed by colliding blobs of the jet 
plasma (e.g., Sikora et al. 1994; Spada et al. 2001; 
Bottcher & Dermer 2010; Mimica & Aloy 2012; Rueda- 
Becerril et al. 2014, see § 5 for the detailed model de¬ 
scription). The model can be however applied to a more 
general situation as well, since it only approximates the 
flaring emission region with a uniform and expanding 
shell of the emitting plasma moving along the outflow, 
and follows self-consistently the evolution of the radiat¬ 
ing particles injected into the shell, but does not involve 
any particular assumption on the exact particle acceler¬ 


ation process involved. 

More specifically, in the model the electron injection 
is assumed to take place at a constant rate with a fixed 
injection spectrum for a given instance of time when the 
shell moves from i?start to i?stop, and to be negligible 
otherwise. The time evolution of the injected electron 
energy distribution is calculated in the emitting plasma 
rest frame (denoted below by primes) following the stan¬ 
dard kinetic equation 


dN^ 

'W 


A 

^7 



T Q7 j 


( 1 ) 


where is the comoving number density of electrons 
with Lorentz factor 7 , and is the electron injection 
function. In this paper we do not prime the electron 
Lorentz factors 7 , noting instead that these always re¬ 
fer to the emitting plasma rest frame. Radiative energy 
losses due to the synchrotron emission and the inverse- 
Compton radiation, as well as the adiabatic losses due to 
the radial expansion of the emitting shell, are all taken 
into account self-consistently in the dy/dt' term. For 
those, the jet magnetic field is assumed to decrease along 
the outflow as B\R) oc corresponding to the case 

of the conserved magnetic energy flux with the dominant 
toroidal component. The energy densities of the external 
target photon fields for the inverse-Compton scattering, 
namely of the BLR and of the hot dusty torus (HDT), 
are assumed to be isotropic in the observer frame and to 
scale with the distance from the core as 

= (2) 

4'7r C1 -|- (Tt/jLextj 

where i?ext is the characteristic radius and Lext repre¬ 
sents the total luminosity of a given field (see Sikora et al. 
2009)^. We note that the BLAZAR code offers the correct 
treatment of the quantum (Klein-Nishina) effects when 


® See in this context Janiak et al. (2015), Appendix Al therein, 
for the alternative scaling of the BLR energy density with the dis¬ 
tance from the central engine. 
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Fig. 1.— Simulated profiles of the analyzed Flare 1 and Flare2 in PKS 1510—089 within the 0.1-300 GeV range (left and right panels, 
respectively), superposed on the 3-hour binned LAT light curve of the source taken from Saito et al. (2013). Simulations denoted in the 
figure by various color curves were performed for different locations of the emitting region along the free-expanding outflow, as explained in 
the text (see §4.1). The values for the corresponding model fits to the decaying phases of the flares are inserted as small figures in the 
panels. The resulting best-fit parameters are given in Table 2. The broad-band emission spectra presented in Figure 3 below coincide with 
the instants marked as “1”, “2”, and “3” in the panels. We note the long-term averaged flux level of the source is 2.7 x 10“^ ph cm”^ s“^ 
above 100 MeV, according to the Fermi-LAT 3FGL catalog (Acero et al. 2015). 


calculating the inverse-Compton emission (see the dis¬ 
cussion in Moderski et al. 2005). 

We parametrize the electron injection spectrum by a 
broken power-law 




1 + 



41 3 (p-g) 


( 3 ) 


source by Nalewajko et al. (2012a). All of these are sum¬ 
marized in Table 2. The jet opening and viewing angles 
were at first assumed as 0jet = dohs = 1/r, where T is the 
bulk Lorentz factor of a moving shell (§4.1); in the fol¬ 
lowing steps of the modeling we also considered smaller 
values of 6jet, in accord with the most recent results of 
the high-resolution radio observations of blazar jets (see 
§4.2 and references therein). 


where Ke is the normalization parameter, p and q are the 
low- and high-energy injection indices, respectively, and 
7 br is the criticial/break electron Lorentz factor. The 
injection function is defined within the range of 7 min < 

T — Tmax- 

After simulating the evolution of the electron energy 
distribution while the shell propagates along the jet 
based on the equation of balance (1), the observed SED 
at a given instance of time is calculated by integrating ra¬ 
diative contributions from all the cells located at different 
radial distances within the shell over the viewing angle 
9ohs- The corresponding 7 -ray light curve is obtained by 
extracting the 7 -ray flux at each step of the simulation 
(see Moderski et al. 2003, for the full description). 

4. TIME-DEPENDENT MODELING 

Several input parameters for the time-dependent sim¬ 
ulations of the analyzed flares were selected based on 
the recent broad-band fitting of the PKS 1510-089 spec¬ 
trum during the 2009 flaring state by Barnacka et al. 
(2014). The main characteristics of the external photon 
fields (BLR and HDT, in particular) were fixed following 
the detailed studies of the accretion disk emission in the 


4.1. “Free-Expanding Jet” Model 

The critical distances along the outflow marking the 
onset and the termination of the electron injection, i?start 
and i?stop, respectively, i.e. the two crucial free pa¬ 
rameters of our modeling, were constrained under the 
condition of a free-expanding outflow Ojet = ^obs = 
1/r = 2.6 deg in the following way. First, the inter¬ 
val Ai? = i?stop — .Rstart was determined from the ob¬ 
served rising time of the flare rg using the general rela¬ 
tion Ai? = cV8Tft/{\-\-z). Next, for a given fixed Ai? we 
varied i?start and i?stop together with the normalization 
of the injection function evaluated the resulting 7 - 
ray fluxes within the 0.1 — 300 GeV range for each set of 
the model parameters (up to the distance i?end > -Rstop), 
and fitted the observed Fermf-LAT light curves with the 
simulated profiles. The final values of the free parame¬ 
ters were then chosen based on the of the model fits 
to the decaying phases of the flares. The results of the 
simulations are presented in Figure 1. 

In the case of the Flare 1, the best model fit to the data 
returns i?start = 0.7 x 10^® cm and i?stop = 0.9 x 10^® cm, 
with the 7 -ray emission settling down around i?end = 
2.3 X 10^® cm with the uncertainty of 0.1 x 10^® cm. Sim- 
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Fig. 2.— The simulated evolution of the electron energy distribution as the emitting shell propagates along the free-expanding jet during 
the analyzed Flare 1 and Flare 2 (left and right panels, respectively), corresponding to the sets of the model parameters given in Table2 . 


ilarly, for the Flare 2 we obtain i?start = 2.3 x 10^® cm, 
l?stop = 3.4 X 10^® cm, and Rend = 6.9 x 10^® cm, with 
the uncertainty of 0.3 x 10^® cm. Overall, the studied 
7 -ray flare light curves can be fitted reasonably well un¬ 
der all the model assumptions specified above, although 
the well-resolved rising profiles of the flares (especially 
of the Flare 2) seem to suggest a more complex time- 
dependence of the electron injection function at the very 
beginning of the acceleration process. 

The simulations performed allow us to propagate the 
evolution of the electron and the broad-band emission 
spectra during the analyzed epochs (though it should 
be noted at the same time that the photon statistics of 
the available Fermi-LAT data precludes us from any pre¬ 
cise characterization of the 7 -ray spectral changes in the 
source in three-hour bins). Figure2 presents the results 
of the simulations regarding the evolution of the electron 
energy distribution as the emitting shell moves along the 
free-expanding outflow. As shown, until the injection 
stops, the number of the injected relativistic electrons 
grows. After the injection terminates, higher energy elec¬ 
trons cool very rapidly, so that the break energy in the 
evolved electron spectrum decreases with time, and the 
high-energy continuum steepens. At the later stages of 
the evolution {R ~ i?end), a moderate pile-up can be 
noted around electron energies 7 < 100. The reason for 
this pile-up is however not the Klein-Nishina suppression 
of the inverse-Compton scattering efficiency (see Dermer 
& Atoyan 2002; Moderski et al. 2005), but instead the 
fact that at such late evolutionary stages the low-energy 
segment of the electron distribution, characterized by the 
injection spectral index p < 2, becomes subjected to the 
efficient radiative cooling dy/dt' oc 7 ^ (see in this con¬ 
text Kardashev 1962). We also note that no pronounced 
spectral hardening due to the Klein-Nishina effects at 
higher electron energies is seen in our simulations, be¬ 
cause the flare is produced effectively outside the BLR 
and the injected high-energy segment of the electron con¬ 
tinuum is steep (c.f. Cerruti et al. 2013); we have how¬ 
ever confirmed that cooled electron spectrum becomes 
softer above electron Lorentz factors 7 ^ 10® when Klein- 


Nishina supression is ignored. 

The corresponding evolution of the broad-band jet 
flaring spectra is shown in Figure 3, including the syn¬ 
chrotron, synchrotron self-Compton (SSC), and the two 
“external-Compton” (EC/BLR and EC/HDT) emission 
components. As shown, in the case of the Flare 1, the 
EC/HDT component dominates the production of the 
high-energy continuum from the observed hard X-ray up 
to ^ 100 MeV photon energies; the soft X-ray range is 
dominated by the SSC emission, the GeV range by the 
EC/BLR process, and in the TeV range the two EC 
components becomes again comparable. In the case of 
the Flare 2, on the other hand, the SSC and EC/BLR 
processes are in general much less relevant, as expected 
taking into account the larger size and distance of the 
emitting region involved in the production of this flare. 

Unfortunately, no simultaneous broad-band data are 
available for the two flares analyzed here. Saito et al. 
(2013) reported that the photon indices within the LAT 
range during the flares were approximately ~ 2 . 0 , which 
is roughly consistent with the spectra simulated in this 
work given the large uncertainty range. For a reference, 
in Figure 3 we included the archival spectral datapoints 
corresponding to the March 2009 flare of PKS 1510-089 
from Barnacka et al. (2014). 

Let us finally comment on the asymmetric profiles of 
the modeled flux enhancements. The standard view on 
this issue is that longer decay phases of blazar flares 
when compared with the flux rising phases reflect longer 
radiative cooling timescales when compared with ex¬ 
tremely rapid, almost instantaneous particle acceleration 
timescales, or even with the short but finite injection in¬ 
tervals controlled by the macroscopic (shock dynamical) 
timescales (as anticipated in our work). However, radia¬ 
tive cooling timescales in PKS 1510-089 jet at sub-pc/pc 
scales are typically expected to be much shorter than 
the observed decay timescales of the 2011 flares (see the 
discussion in Saito et al. 2013). 

And indeed, our simulations reveal that the dominant 
factor shaping the flux decay profiles is not the radia¬ 
tive cooling of the high-energy particles, but instead a 




6 


Saito et al. 




Fig. 3.— The simulated time evolution of the broad-band SEDs of PKS 1510-089 during the analyzed Flare 1 and Flare2 (left and right 
panels, respectively). The spectra were extracted at the instants marked as “1”, “2”, and “3” in Figure 1 (black, red, and green curves, 
respectively). Vertical dashed lines indicate the energy range of 100 MeV - 300 GeV. For a reference, in the plots we include also the 
archival spectral datapoints corresponding to the March 2009 flare of the source from Barnacka et al. (2014). Calendar dates are October 
19.5-20.5, 2011 for Flare 1 and November 7-9, 2011 for Flare 2. 

gradient of the Doppler factor across the emitting shells 
of the jet plasma. Namely, for the model parameters 
considered above in this section, the emitting shells are 
relatively thin but instead considerably extended in the 
radial direction (see §5), so that different parts of the 
shells are observed at different viewing angles. As a re¬ 
sult, emission produced within the parts characterized by 
the largest inclinations arrive to the observer with a sig- 
nihcant delay when compared with the emission output 
of the parts located at smaller viewing angles. 

This effect is visualized in Figured, where we present 
the simulated profiles of the Flare 1 corresponding to the 
different values of the jet opening angle (equal by as¬ 
sumption to the jet viewing angle), and all the other 
model parameters fixed as before except for the adjusted 
electron normalization. As shown, with the decreas¬ 
ing jet opening angle the variance of the Doppler factor 
across the emitting shell becomes smaller, and as a result 
the flare’s asymmetry decreases. Only for the dramati¬ 
cally small djet = ^obs = 0.3 deg the simulated decaying 
timescale becomes comparable to the radiative cooling 
timescale of the 7 -ray emitting electrons. 

4.2. “Collimated Jet” Model 

Due to the causality requirement, for a free-expanding 
jet one has F 0 jet — 1 at most, the condition which is 
often anticipated in blazar modeling. However, several 
recent radio studies of relativistic jets in blazar sources 
imply highly collimated outflows on milli-arcsec scales, 
with small opening angles F 0jet ~ 0.1 (e.g., Clausen- 
Brown et al. 2013; Jorstad et al. 2005; Zdziarski et al. 

2015). In order to investigate such a possibility in more 
detail, we repeated simulations of the selected flares in 
the source assuming different values of 0.1 < F 0jet < 1 
(for the fixed F = 22). We again imposed the condi¬ 
tion 0 obs = ^jet, so that the relativistic beaming is max¬ 
imized; this requirement is necessary, since in the case 


of a misaligned PKS 1510-089 jet (i.e., for 0obs > ^jet), 
unrealistically large electron injection is needed in order 
to produce the observed 7 -ray flux enhancements. 

The resulting best-fit positions of the onset of the 7 - 
ray emitting regions for different values of ^jet) and the 
corresponding normalizations of the electron injection, 
are presented in Figure 5. For both flares considered, the 
location of the emission zone increases further away from 
the SMBH as the jet opening angle decreases. This de- 
pendance may be understood by considering the effect 
of the Doppler factor gradient discussed in the previ¬ 
ous section. In particular, for smaller and smaller 0 jet, 
the Doppler factor variance across the emitting shell de¬ 
creases; as a result, the emission region have to be placed 
further and further away from the jet base so that the 
increased electron cooling timescales can support the ex¬ 
tended decay phases of the flares (note that at larger 
distances from the core the energy densities of the BLR 
and HDT photon fields decrease, and hence the corre¬ 
sponding radiative cooling timescales increase). 

In the case of the highly collimated jet with particu¬ 
larly small 0]et — 0 .26 deg, the best-fit position of the 
emission zone turns out as i?start — 4.3 x 10^® cm for the 
Flare 1, and i?start — 8.1 x 10^® cm for the Flare 2. To¬ 
gether with the simulation results presented in the pre¬ 
vious section, this indicates therefore a relatively nar¬ 
row range allowed for the production of the observed 
7 -ray flares in PKS 1510-089, namely 10^® cm > Rstart > 
10^® cm for the jet opening angle 0.1 < F 0jet < 1- At the 
same time, the required amount of the electron injection 
for producing the observed 7 -ray luminosity may be sub¬ 
stantially smaller for a well-collimated outflow (see the 
lower panels in Figure 5). This means that in the case of 
a jet with F 0jet < 1 it is “easier” to produce huge 7 -ray 
outbursts with the total released radiative power of the 
order of the accretion disk luminosity, while somewhat 
extreme conditions are needed for such a maximum effi- 

























7 


Modeling of Gamma-ray Flares in PKS 1510-089 


Bjet = 2.6° Bjet = 0.26° 





Fig. 4.— Simulated profiles of the Flare 1 assuming different val¬ 
ues of the jet opening angle (and other parameters as given in 
Table2 except for the adjusted electron normalization), which il¬ 
lustrate the effect of the Doppler factor gradient across the emitting 
shell in the PKS 1510-089 jet (see §4.1 for the discussion). 


Fig. 5.— The best-fit location of the onset of the flaring zone 
(upper panels) and the corresponding normalization of the electron 
injection function (lower panels) for the Flare 1 and Flare 2 (left and 
right panels, respectively), as functions of the jet opening angle 
(see §4.2 for the discussion). 


ciency in the case of a standard, free-expanding jet with 
rdjet — 1 (see the discussion in the following section). 


5. DISCUSSION 

In the framework of the internal shock scenario, a col¬ 
lision of two plasma blobs characterized by different bulk 
Lorentz factors Fi = (1 — and F 2 = (1 — 

produces the double-shock shock structure propagating 
within the merged portion of the jet matter. It is as¬ 
sumed that the electrons are accelerated at the shock 
fronts, and are injected with a given energy spectrum 
into the downstream region where they cool radiatively 
and adiabatically. The resulting shock velocity can be 
found under the specific assumptions on the dynamics 
of the colliding plasmoids. In particular, if the collid¬ 
ing portions of the jet flow differ only in bulk veloc¬ 
ities, and the jet magnetic field is negligible dynami¬ 
cally, as assumed in this paper, the downstream (emit¬ 
ting) plasma is characterized by the bulk Lorentz factor 

F = (1 — ~ VFi 62 as long as F 2 > Fi 3> 1 

(e.g., Stawarz et al. 2004; Moderski et al. 2004). In the 
upstream region rest frame (denoted below by double¬ 
primes), moving with the velocity /3i relative to the ob¬ 
server, the bulk Lorentz factor of the shock is 


where 


F" = 

^ sh 


13" = 


-kl)(4F"-l)2 

( 4 ) 

8 F" -k 10 

13-Pi 
l-/3/3i ■ 

( 5 ) 


This implies that the resulting shock velocity in the 
downstream region rest frame (denoted in this paper by 


primes), 

saturates at the approximately constant value of /3'j, ~ 
0.1 (Stawarz et al. 2004). 

The observed rising timescale of a flare, ta, may be 
identified with the time interval rinj when active shocks 
propagating through the outflow inject freshly acceler¬ 
ated electrons into the downstream region (which is a 
moving portion of the jet matter). Hence, the shell width 
in the emitting region rest frame is 

C = 2c/3'hT-4j:^2c/i:h<5r9. (7) 

On the other hand, the transverse extension of the emit¬ 
ting region can be estimated simply as 

r' = r ~ Rstop Ojet ■ ( 8 ) 

In the case of the Flare 1, for which the observed rising 
timescale is of the order of three hours, one has for exam¬ 
ple ~ 1.5x 10^^ cm and r' ~ 4x 10^® cm under the con¬ 
dition of a free-expanding jet with djet = dobs = 1/F ~ 
2.6 deg; the emerging ratio ^ 0.035 implies there¬ 

fore an extremely thin shell of the emitting plasma, even 
at the final stages of the shock evolution. Note, however, 
that in the observer frame an almost spherical apparent 
shape would be measured, namely £sh/r = S£'^^/r' ~ 0 . 8 . 
The intrinsic aspect ratio of the emitting region increases, 
on the other hand, in the case of a well-collimated jet. 
In particular, with the jet opening and viewing angles 
djet ~ dobs — 0.26 deg, the downstream region linear scale 
reads as ~ 3 x 10^® cm, and the emitting region ra¬ 
dius becomes r' ~ 4 x 10^® cm, leading to £si^/r' ^ 0.7, 
or £sh/i’ ~ 30 in the observer frame. In the case of 
the Flare 2, characterized by the observed flux rising 
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timescale of the order of 12 hours, the analogous val¬ 
ues read as ^ 6 x 10^® cm and r' ~ 1.5 x 10^^ cm for 
0jet — 2.6 deg, or ~ 10^® cm and r' ~ 1.5 x 10^® cm 
for 0jet — 0.26 deg, so that the corresponding aspect ratio 
are basically the same as for the Flare 1. 

The total number of electrons injected into the emit¬ 
ting shell during the shock activity period, 

M; = y it' J dy = T-'^j X JdjQ^, (9) 


together with the electron mean energy, 


(7inj) 


fdjlQy 

fdj 


( 10 ) 


gives the jet comoving energy density of the 
electrons 

, _ (binj) 


radiating 


( 11 ) 


Hence, the electron kinetic flux (including also the power 
dissipated radiatively during the flare) reads as 


Le = 7rr^ c/3r^ u'^ ( 12 ) 


(see equation 7). With the best-fit parameters of the free- 
expanding jet model (see Table 2), one obtains therefore 
Le — 1.3 X lO'^’^ergs”^ for the Flare 1, and Le — 0.5 x 
10^^ergs“^ for the Flare 2. 

With the aforementioned model parameters, the elec¬ 
tron mean energy reads as (yinj) — HO. In the frame¬ 
work of the anticipated internal shock scenario, assum¬ 
ing in addition that the upstream jet plasma is cold 
(meaning at most trans-relativistic energies of plasma 
particles), this value excludes a pure electron-positron 
content of the PKS 1510-089 jet. That is because in 
the case of a pure pair plasma, the energy conservation 
(binj) rWe TTieC^ ~ iF 2 WemgC^ would then imply un¬ 
realistically high bulk Lorentz factor of the faster shell, 
namely r 2 ~ 5000 (for the the downstream bulk Lorentz 
factor F = 22). Hence, self-consistency of the modeling 
presented in this paper requires protons to dominate the 
plasma inertia. 

Let us therefore assume that the colliding shells with 
total energies Ei and E 2 are initially dominated by 
cold protons with the total number Afp. Under the 
r 2 > Fi ^ 1 assumption specified previously, one can 
derive the efficiency of the energy dissipation after the 
shells’ collision as 


^diss — 


Ediss (r2/r - if 


El + E2 (r2/r)^ -h 1 


(13) 


Quantifying next the energy dissipated per proton in the 
TTipC^ units as 


K = 


E'diss 

Afp nipC^ ’ 


(14) 


and keeping in mind that Ei = ^TiAfpiripC^, E 2 = 
^T 2 AfpmpC^, and = Udiss/P, one can find 


^ 1 (r 2 /r - 1 )^ 
“ 2 r 2 /F 


(15) 


And since {'jinf Afe meC^ = ?7e A'diss^ where rje is the effi¬ 
ciency of the energy transfer to relativistic electrons at 
the shock front, the composition of the jet can be finally 
estimated as 


Afe mp me 

TT - "He K, ^ 

ffv (7inj) 


(16) 


(see, e.g., Moderski et al. 2004). The typically consid¬ 
ered value of ? 7 e — 0.5 and a pure electron-proton jet 
composition Afg/Afp ~ 1 imply therefore r 2 ~ 35 with 
the corresponding ?7diss — 0.1. Any larger amount of 
electron-positron pairs in the jet would increase both the 
bulk Lorentz factor of the faster shell (for the given F), 
and the overall energy dissipation efficiency; for example, 
Afe/Afp ~ 3 (meaning the pair content Afe+IAfe- — 0-5) 
would require a still reasonable value of r 2 ~ 50, with 
the corresponding r/diss — 0.25. Note in this context that 
the outflows remains dominated dynamically by protons. 


Lp 1 mp/me 

Le 'Afe/Afp (^inj) 


(17) 


as long as Me/Afp < 10. In particular, with Afe/Afp — 
1 — 3 adopted hereafter, one has Lp/Le ^ 10. 

The total luminosity of the accretion disk in PKS 1510- 
089 was estimated by Nalewajko et al. (2012a) as Ldisk — 
5 X 10"^^ ergs“^, meaning the accretion power in the sys¬ 
tem of the order of Lace ~ lO^^ergs”^ (for the assumed 
the standard, 10% radiative efficiency factor). The model 
parameters evaluated above under the F 6-^et — 1 con¬ 
dition, imply therefore an extreme efficiency of the jet 
production, with the total jet kinetic power Lj Lp 
10 Lace- In the case of a highly collimated jet, this ef- 
hciency may be decreased quite significantly (see Fig¬ 
ure 5). On the other hand, the jet magnetic held within 
the blazar emission zone, which is rather weak already 
in the free-expanding jet case, becomes even less rele¬ 
vant dynamically with the decreasing product FSjet- In 
particular, with 0jet — 2 .6 deg, the ratio of the electron 
kinetic energy and Poynting huxes. 


Le 2meC^ fdjjQj 

^ " c/3(h "" lR0jetB'(R)f ’ 


reads as Le/LB ~ 60 and ^ 20 for the Flare 1 and the 
Flare 2, respectively. Assuming instead 0jet — 0.26 deg, 
but keeping a comparably large Compton dominance in 
the source (i.e., the ratio of the high-energy and syn¬ 
chrotron peak luminosities ~ 10 — 100), we obtain wor- 
risomely small jet magnetization of Le/LB ^ 100. 

It is interesting to note in this context that the 
model ht to the average/quiescence broad-band spec¬ 
trum of PKS 1510-089 using the BLAZAR code, presented 
in Kataoka et al. (2008), indicates Lp ~ 2 x 10^®ergs“^ 
and Le — Lb — 0.6 x 10^®ergs“^. This, together with 
the results of our simulations, may suggest that while 
during the extended source quiescence the outflow is in 
equipartition between the relativistic electron and mag¬ 
netic field energy fluxes, with the total jet kinetic lu¬ 
minosity dominated only slightly by cold protons and 
constituting only some smaller fraction of the available 
accretion power (Lace ^ 10Lp with Lg ^ Lb ~ 0.3Lp), 
dramatic though relatively rare flaring events consist of 
an excess energy flux carried out predominantly by the 
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Gjet = 2.6° Gjet = 0.26° 



Fig. 6. — Simulated light curves corresponding to the Flare 1 in 
optical, X-rays, HE 7 -rays, and VHE 7 -rays (top to bottom, respec¬ 
tively), for the cases of a free expanding conical jet ( 0 jet = ^obs = 
1/F, left panels), and a collimated jet (with ^jet — ^obs — 0 - 1 /r, 
right panels). The location of the emission zone and the normal¬ 
ization of the electron injection are fixed at the best-fit values for 
the HE 7 -ray flares (see Figure 4). Various emission components 
are denoted by different styles of the curves (dashed, dot-dashed, 
dot-dot-dashed, and dotted for the synchrotron, SSC, EC/BLR, 
and EC/HDT, respectively), with the solid curves corresponding 
to the sum of all the emission components in each panel. 

jet particles, and exceeding the mean accretion power in 
the source {Lp ~ lOLacc and Lp ~ lOig ^ Lb)- 
The above considerations should be taken with ex¬ 
treme caution, however, since in the modeling presented 
here no synchrotron data concurrent with the analyzed 
7 -ray flares could be utilized. Any robust estimation of 
the magnetic field within the emitting region would, in 
fact, require exactly simultaneous infrared/optical and 
GeV flux measurements. At the same time, we note that 


Gjet = 2.6° Gjet = 0.26° 



Fig. 7. — Same as Figure 6, but for Flare 2. 


small jet opening angle F 0jet < 1 implies higher number 
density of the radiating electrons when compared with 
the case of a free-expanding jet, regardless on the jet 
magnetization, and therefore an elevated SSC spectral 
component (which may dominate the entire X-ray do¬ 
main, up to even the soft y-ray regime). Hence, high- 
quality X-ray data simultaneous with the GeV flaring 
events, when modeled as presented in this paper, could, 
in principle, constrain robustly the jet opening angle in 
the source. This is demonstrated more quantitatively in 
Figures 6 and 7, where we present the simulated light 
curves corresponding to the Flare 1 and Flare 2, respec¬ 
tively, in various frequency ranges, for different jet open¬ 
ing angles. One can see that flaring timescales become 
in general shorter in the case of a collimated outflow, the 
most pronouncedly however in the X-ray domain, which 
constitutes a complex superposition of the SSC and EC 
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emission components. Also, for the given HE 7 -ray flare 
profile and amplitude there is a significant difference be¬ 
tween the relative amplitudes of optical/X-ray and VHE 
7 -ray flares depending on the jet opening angle; in par¬ 
ticular, in the case of ^jet = 1/r the optical and X-ray 
flares are much more prominent when compared with a 
moderate increase in the VHE flux, while for Oiet = 0 . 1 /r 
the situation reverse. 

6 . CONCLUSIONS 

In this paper we presented a new approach for con¬ 
straining luminous blazars, incorporating fully time- 
dependent and self-consistent modeling of particularly 
bright and well-resolved 7 -ray flares of PKS 1510-089 de¬ 
tected with Fermi-LAT. The two flares selected for the 
analysis, studied before by Saito et al. (2013), consti¬ 
tute the best known examples of prominent, isolated, 
and coherent events (with well-defined flux rising and de¬ 
cay phases on hourly timescales), unlike the majority of 
the observed blazar flux enhancements which seem rather 
like a superposition of distinct (though possibly related) 
but just unresolved sub-flaring/flickering. Unfortunately, 
no simultaneous data at lower (radio-to-X-ray) frequen¬ 
cies are available for the analyzed flares. 

The results of our modeling, performed with the 
BLAZAR code developed by Moderski et al. (2003) in the 
framework of the internal shock scenario, are largely in 
agreement with the recent broad-band observations of 
FSRQs in general, and with the detection of TeV 7 - 
ray photons correlated with the GeV flares in particu¬ 
lar. Such a correlation, along with the apparent smooth¬ 
ness of the observed 7 -ray spectra from 100 MeV up to 
the TeV range, strongly suggests a co-spatiality of the 
GeV and TeV emitting regions, which have to be in ad¬ 
dition located outside the BLR in order to avoid a sig- 
nihcant attenuation of the 7 -ray fluxes due to the ef¬ 
ficient photon-photon annihilation (see Barnacka et al. 
2014). And indeed, the best-fit location of the 7 -ray flar¬ 
ing region estimated here for PKS 1510-089 turns out as 
R ~ 0.3 pc for a free-expanding jet with the opening an¬ 
gle 6»jet ~ l/P ~ 2.6deg, up to i? ~ 3pc for a collimated 
outflow with 0 jet — 0 .l/P ~ 0.26deg (as advocated by 
Clausen-Brown et al. 2013; Jorstad et al. 2005; Zdziarski 
et al. 2015). This is safely beyond the characteristic scale 
of the BLR in the source (^ 0.03pc). 

We note that several other complementary arguments 
have been presented in the literature in support of the 
dominant blazar emission zone located at ~ pc distances 
from the core (see the discussion in Sikora et al. 2009; 
Nalewajko et al. 2014). 

Under the P 0jet — 1 assumption, our modeling in¬ 


Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, ApJ, 699, 
817 

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2011, ApJ, 733, 
L26 

Acero, F., Ackermann, M., Ajello, M., et al. 2015, 
arXiv:1501.02003v2 

Abramowski, A., Acero, F., et al. 2013, A&A, 554, AA107 
Aleksic, J., Ansoldi, S., Antonelli, L. A., et al. 2014, A&A, 569, 
A46 

Asano, K., Takahara, F., Kusunose, M., Toma, K., & Kakuwa, J. 
2014, ApJ, 780, 64 


dicates in addition an extremely efficient jet produc¬ 
tion during the flaring events and a mixed jet content, 
with the dominant proton energy flux exceeding the 
total available accretion power in the studied blazar, 
Lj Lp 10 Lace- This is in contrast to the quies¬ 
cence states of the source, during which the Lj ~ 0.1 Lace 
condition and an approximate equipartition between dif¬ 
ferent plasma constituents (protons, electrons, and mag¬ 
netic field) seem to hold. In the case of a collimated 
jet with P 0 jet — 0 . 1 , on the other hand, the flaring 
jet production efficiency decreases by an order of mag¬ 
nitude. Only strictly simultaneous observations of flar¬ 
ing PKS 1510-089 at infrared. X-ray, and GeV photon 
energies, augmented by fully self-consistent and time- 
dependent simulations as presented in this paper, may 
help to remove such an order-of-magnitude uncertainty, 
by imposing precise constraints on the magnetization and 
opening angle of the emitting region. 

We note that the jet production efficiency Lj/Lacc ex¬ 
ceeding 100 % is in principle consistent with the recent 
understanding of the jet launching via the Blandford- 
Znajek process in the case when the magnetic field 
threading the SMBH horizon saturates at the maximum 
sustainable level (see Tchekhovskoy et al. 2011; McKin¬ 
ney et al. 2012, also the discussion in Ghisellini et al. 
2014). Yet the jet magnetization emerging in our model¬ 
ing is very low, so a self-consistency of the internal shock 
scenario explored in this paper would require in addition 
a very efficient conversion of the magnetic energy flux 
to the particle flux (accompanied by an efficient proton 
loading) between the jet base and the blazar emission 
zone (located at ~pc distances from the core). It may 
be that other energy dissipation processes are involved 
instead in the production of blazar flares, including tur¬ 
bulent acceleration (e.g., Ushio et al. 2010; Lefa et al. 
2011; Tramacere et al. 2011; Cao & Wang 2013; Yan et 
al. 2013; Asano et al. 2014; Zheng et al. 2014; Kakuwa 
et al. 2015; Chen et al. 2015), or relativistic magnetic 
reconnection (see, e.g., Narayan & Piran 2012; Giannios 
2013; Sironi et al. 2015, and references therein). These 
interesting alternatives have been analyzed so far more 
quantitatively only in the context of low-power blazars 
of the BL Lac type, and not FSRQs, for which the inter¬ 
nal shock scenario is still the most appealing possibility 
(though see in this context also the discussion in Nale¬ 
wajko et al. 2012 b). 
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